#!/usr/bin/perl
use strict;

#my ($RANGE, $CLUSTERS, $MIN_CLUSTER_SIZE, $ADDITIONAL_CNT) = (300, 30, 10, 300);
my ($RANGE, $CLUSTERS, $MIN_CLUSTER_SIZE, $ADDITIONAL_CNT) = (100, 3, 10, 20);
my $JUMBLE = 0;

sub sort_arr_of_arrays
{
	for my $sub_arr(@_)
	{
		@$sub_arr = sort @$sub_arr;
	}
	@_ = sort {$a->[0] <=> $b->[0] } @_;
	
}


#
# Divide indices into clusters
#
print STDERR "Divide indices into clusters\n";
my @clusters;
while (1)
{
	for (0.. ($RANGE - 1))
	{
		push(@{$clusters[int(rand($CLUSTERS))]}, $_);
	}

	# make sure no cluster is too small
	for (@clusters)
	{
		if (@$_ < $MIN_CLUSTER_SIZE)
		{
			@clusters = ();
			redo;
		}
	}
	last;
}

#
# Create pairs joining up clusters
#
print STDERR "Create pairs joining up clusters\n";
my @pairs;

#pairs saved by cluster
my @pairs_by_cluster_mine;

for my $cluster_id(0.. ($CLUSTERS - 1))
{
	my $cluster = $clusters[$cluster_id];
	#
	# pairs should make up 75% of possible pairs for this cluster
	#
	my $clust_size = @$cluster;
	my $max_pairs = int($clust_size * ($clust_size - 1) * 3 / 8);
	
	# matrix to prevent redundant pairs
	my %paired;
	
	for my $pair_cnt(0..$max_pairs)
	{
		my $id1 = $cluster->[int(rand($clust_size))];
		my $id2 = $cluster->[int(rand($clust_size))];

		# no self-self or redundant pairs
		redo if $id1 == $id2;
		redo if (exists $paired{$id1}{$id2});
		$paired{$id1}{$id2} = 
		$paired{$id2}{$id1} = 1;

		# first position is random number for jumbling up 
		my $pair_str = "$id1\t$id2\t0\n";
		push(@pairs,  [rand(), $pair_str]);
		push(@{$pairs_by_cluster_mine[$cluster_id]}, $pair_str);
	}
}

#
#	jumble pairs
#
if ($JUMBLE)
{
	@pairs = sort {$a->[0] <=> $b->[0]} @pairs;
}


#
#	print pairs
#
print STDERR "Print pairs\n";
{
	open TEST, ">original.pairs" or die;
	print TEST "seq1\tseq2\n----\t----\n";
	print TEST $_->[1] for (@pairs);
}

#
#	Read clustered pairs back in 
#
print STDERR "Read clustered pairs back in\n";
my @pairs_by_cluster_bin;
{
	my @clustered_pairs =`seq_pairs_cluster_by_single_linkage -v --input_file original.pairs`;
	system("seq_pairs_cluster_by_single_linkage -I --input_file original.pairs > original.identifiers");
	my $cnt_clusters_bin;
	for (@clustered_pairs)
	{
		if (/^>/)
		{
			$cnt_clusters_bin++;
			next;
		}
		push(@{$pairs_by_cluster_bin[$cnt_clusters_bin - 1]}, $_);
	}
}

#
#	sort and compare clustered pairs with original clustered pairs
#
print STDERR "Sort and compare clustered pairs with original clustered pairs\n";
{
	sort_arr_of_arrays(@pairs_by_cluster_bin);
	sort_arr_of_arrays(@pairs_by_cluster_mine);
	
	if (@pairs_by_cluster_mine == @pairs_by_cluster_bin)
	{
		print STDERR "Success\n\n";
	}
	else
	{
		print STDERR "Error\n\n";
	}
}















#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#
#
#	Create additional items
#


sub scramble
{
	return map {$_->[1]}					# remove first entry
				sort{$a->[0] <=> $b->[0]}	# sort by random number
				map {[rand(), $_]}			# make first entry random number
				@_;
}

sub add_pairs_to_cluster($\@$\@$)
{
	my ($cluster, $extra_ids, $pairs_by_cluster_mine, $pairs, $value) = @_;

	#
	#	Create pairs
	#
	# scramble order of extra_ids
	my @cluster_members = scramble(@$cluster);

	#throw away half
	splice(@cluster_members, @cluster_members / 2);

	my $cnt_additional_pairs = 0;
	# add pairs for rest
	for my $id1(@$extra_ids)
	{
		for my $id2(@cluster_members)
		{
			# first position is random number for jumbling up 
			my $pair_str = "$id1\t$id2\t$value\n";
			push(@$pairs,  [rand(), $pair_str]);
			push(@$pairs_by_cluster_mine, $pair_str);
			++$cnt_additional_pairs;
		}
	}
}

my $additional_identifiers = 0;
my @parachuted_clusters;
my $curr_id = $RANGE;
for my $cluster_id(0.. ($CLUSTERS - 1))
{
	my $cluster = $clusters[$cluster_id];

	#
	#	Create extra $ADDITIONAL_CNT items per cluster
	#
	my @extra_ids;
	push(@extra_ids, $curr_id++) for (0.. ($ADDITIONAL_CNT - 1));
	$additional_identifiers += $ADDITIONAL_CNT;

	# save extra ids
	push(@{$parachuted_clusters[$cluster_id]}, @extra_ids);
	
	#
	#	Create pairs
	#
	add_pairs_to_cluster($cluster, @extra_ids, $pairs_by_cluster_mine[$cluster_id], 
						 @pairs, 1);

}

#
#	Create items in both cluster #0 and cluster #1
#
{
	die unless $CLUSTERS > 1;
	#
	#	Create extra $ADDITIONAL_CNT items per cluster
	#
	my @extra_ids;
	push(@extra_ids, $curr_id++) for (0..($ADDITIONAL_CNT - 1));
	$additional_identifiers += $ADDITIONAL_CNT;

	# save extra ids
	push(@{$parachuted_clusters[0]}, @extra_ids);
	push(@{$parachuted_clusters[1]}, @extra_ids);

	#
	#	Create pairs
	#
	add_pairs_to_cluster($clusters[0], @extra_ids, $pairs_by_cluster_mine[0], 
						 @pairs,2);
	add_pairs_to_cluster($clusters[1], @extra_ids, $pairs_by_cluster_mine[1], 
						 @pairs,2);
}


#
#	Create X-link between cluster #0 and cluster #1
#		these are discarded
#
{
	die unless $CLUSTERS > 1;
	my $cnt_cross_links = 0;
	my $cluster0 = $clusters[0];
	my $cluster1 = $clusters[1];
	for my $id1(@$cluster0)
	{
		for my $id2(@$cluster1)
		{
			# first position is random number for jumbling up 
			my $pair_str = "$id1\t$id2\t-2\n";
			push(@pairs,  [rand(), $pair_str]);
			last if (++$cnt_cross_links > $ADDITIONAL_CNT);
		}
	}
	print STDERR "\t", ($cnt_cross_links), 
				" pairs should be ignored because they x-link pre-existing clusters.\n";
}



#
#	Create items in both cluster #0 and cluster #1 && which are linked to one another
#
my @extra_ids;
{
	die unless $CLUSTERS > 1;
	#
	#	Create extra $ADDITIONAL_CNT items
	#
	push(@extra_ids, $curr_id++) for (0..($ADDITIONAL_CNT - 1));
	$additional_identifiers += $ADDITIONAL_CNT;

	# save extra ids
	push(@{$parachuted_clusters[0]}, @extra_ids);
	push(@{$parachuted_clusters[1]}, @extra_ids);
	
	#
	#	Create pairs
	#
	add_pairs_to_cluster($clusters[0], @extra_ids, $pairs_by_cluster_mine[0], 
						 @pairs,4);
	add_pairs_to_cluster($clusters[1], @extra_ids, $pairs_by_cluster_mine[1], 
						 @pairs,4);

	# add pairs for links between parachuted items
	# these should be in both clusters
	my $cnt_in_multiple_clusters = 0;
	for my $id1(@extra_ids)
	{
		for my $id2(@extra_ids)
		{
			next if $id1 == $id2;
			# first position is random number for jumbling up 
			my $pair_str = "$id1\t$id2\t8\n";
			push(@pairs,  [rand(), $pair_str]);
			push(@{$pairs_by_cluster_mine[0]}, $pair_str);
			push(@{$pairs_by_cluster_mine[1]}, $pair_str);
			++$cnt_in_multiple_clusters;
		}
	}
	print STDERR "\t", $cnt_in_multiple_clusters * 2,
		" pairs should be added because they are present in 2 clusters.\n";
}

#
#	Create items in cluster #0 and cluster #1 && which are linked to one another
#
{
	die unless $CLUSTERS > 1;
	#
	#	Create extra $ADDITIONAL_CNT items
	#
	my @extra_ids1;
	push(@extra_ids1, $curr_id++) for (0..($ADDITIONAL_CNT - 1));
	$additional_identifiers += $ADDITIONAL_CNT;
	my @extra_ids2;
	push(@extra_ids2, $curr_id++) for (0..($ADDITIONAL_CNT - 1));
	$additional_identifiers += $ADDITIONAL_CNT;

	# save extra ids
	push(@{$parachuted_clusters[0]}, @extra_ids1);
	push(@{$parachuted_clusters[1]}, @extra_ids2);
	
	#
	#	Create pairs
	#
	add_pairs_to_cluster($clusters[0], @extra_ids1, $pairs_by_cluster_mine[0], 
						 @pairs,4);
	add_pairs_to_cluster($clusters[1], @extra_ids2, $pairs_by_cluster_mine[1], 
						 @pairs,4);

	# add pairs for links between parachuted items
	# these should be ignored
	my $cnt_in_multiple_clusters = 0;
	for my $id1(@extra_ids1)
	{
		for my $id2(@extra_ids2)
		{
			# first position is random number for jumbling up 
			my $pair_str = "$id1\t$id2\t8\n";
			push(@pairs,  [rand(), $pair_str]);
			++$cnt_in_multiple_clusters;
		}
	}
	print STDERR "\t", $cnt_in_multiple_clusters,
		" pairs should be ignored because they are parachuted into 2 different clusters.\n";
}

print STDERR "\t", $additional_identifiers,
				" identifiers parachuted into pre-existing clusters.\n";


#
#	Spurious self-self orphans
#
{
	#
	#	Create extra $ADDITIONAL_CNT items
	#
	my @self_self_orphans;
	push(@self_self_orphans, $curr_id++) for (0.. ($ADDITIONAL_CNT - 1));

	for my $id1(@self_self_orphans)
	{
		# first position is random number for jumbling up 
		my $pair_str = "$id1\t$id1\t-3\n";
		push(@pairs,  [rand(), $pair_str]);
	}
	print STDERR "\t", $ADDITIONAL_CNT,
				" self - self pairs should be ignored.\n";
}



	
#
#	Create extra items altogether which are linked to new cluster only
#
my @orphan_ids;
{
	#
	#	Create extra $ADDITIONAL_CNT items
	#
	push(@orphan_ids, $curr_id++) for (0.. ($ADDITIONAL_CNT - 1));

	
	# add pairs for links between parachuted items
	# these should not be in any preexisting cluster
	for my $id1(@orphan_ids)
	{
		for my $id2(@extra_ids)
		{
			# first position is random number for jumbling up 
			my $pair_str = "$id1\t$id2\t-1\n";
			push(@pairs,  [rand(), $pair_str]);
		}
	}

	my $cluster_id = @pairs_by_cluster_mine;

	# add pairs for orphans
	# these should be in a new cluster
	for my $id1(@orphan_ids)
	{
		for my $id2(@orphan_ids)
		{
			next if $id1 == $id2;

			# first position is random number for jumbling up 
			my $pair_str = "$id1\t$id2\t-1\n";
			push(@pairs,  [rand(), $pair_str]);
			push(@{$pairs_by_cluster_mine[$cluster_id]}, $pair_str);
		}
	}
	print STDERR "\t", ($ADDITIONAL_CNT * $ADDITIONAL_CNT), 
				" pairs should be ignored because only one has been parachuted\n";
}

#
#	jumble pairs & print
#
print STDERR "Print extra pairs\n";
{
	if ($JUMBLE)
	{
		@pairs = sort {$a->[0] <=> $b->[0]} @pairs;
	}
	open TEST, ">extra.pairs" or die;
	print TEST $_->[1] for (@pairs);
}

#
#	Write out extra ids
#
{
	open(EXTRA_IDS, ">extra.identifiers") or die;
	for my $cluster_id(0.. ($CLUSTERS - 1))
	{
		my $ids			= $clusters[$cluster_id];
		my $extra_ids	= $parachuted_clusters[$cluster_id];
		my @total_ids	= @$ids;
		push(@total_ids, @$extra_ids) if ($extra_ids && @$extra_ids);
		print EXTRA_IDS (join ("\t", sort (@total_ids)), "\n");
		print EXTRA_IDS ("\tOriginal:\t", join ("\t", sort @$ids), "\n");
		print EXTRA_IDS ("\tExtra:\t", join ("\t", sort @$extra_ids), "\n") if ($extra_ids && @$extra_ids);
	}
	print EXTRA_IDS ("Orphaned IDs:\n", join ("\t", sort @orphan_ids), "\n");
}


#
#	Read clustered pairs back in 
#
print STDERR "Cluster and parachute using seq_pairs_cluster_by_single_linkage...\n";
my @pairs_by_cluster_bin;
{
	my @clustered_pairs =`seq_pairs_cluster_by_single_linkage -v -p original.identifiers < extra.pairs`;
	system("seq_pairs_cluster_by_single_linkage -vvv -p original.identifiers < extra.pairs > clustered.pairs");
	system("seq_pairs_cluster_by_single_linkage -I -p original.identifiers < extra.pairs > clustered.identifiers");
	system("seq_pairs_cluster_by_single_linkage -I < extra.pairs > clustered.without_parachuting.identifiers");
	my $cnt_clusters_bin;
	for (@clustered_pairs)
	{
		if (/^>/)
		{
			$cnt_clusters_bin++;
			next;
		}
		push(@{$pairs_by_cluster_bin[$cnt_clusters_bin - 1]}, $_);
	}
}

#
#	sort and compare clustered pairs with original clustered pairs
#
print STDERR "Sort and compare clustered pairs with original clustered pairs\n";
{
	sort_arr_of_arrays(@pairs_by_cluster_bin);
	sort_arr_of_arrays(@pairs_by_cluster_mine);
	
	if (@pairs_by_cluster_mine == @pairs_by_cluster_bin)
	{
		print STDERR ('*' x 60), "\n\n";
		print STDERR "Success\n\n";
		unlink(glob("*.pairs"));
		unlink(glob("*.identifiers"));
		print STDERR ('*' x 60), "\n\n";
	}
	else
	{
		print STDERR ('!' x 60), "\n\n";
		print STDERR ('!' x 60), "\n\n";
		print STDERR ('!' x 60), "\n\n";
		print STDERR "Error\n\n";
		print STDERR ('!' x 60), "\n\n";
		print STDERR ('!' x 60), "\n\n";
		print STDERR ('!' x 60), "\n\n";
		print STDERR ('!' x 60), "\n\n";
	}
}




